//////////////////////////////////////////////////////////////////////////
/// Copyright (c) 2011, Sidorenko 'PC' Alexander                       ///
/// mailto: sidorenko /*dot*/ alexander (at) gmail [dot] com           ///
/// This file is a part of the Radiating Gas Dynamics project          ///
//////////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "SodProblem.h"
#include "ColorBar.h"

#include "PlotHelper.h"

const gdScalar T0 = gdScalar(2273.0);
const gdScalar T1 = gdScalar(293.0);

const gdScalar p0 = gdScalar(1.01e5);
const gdScalar p1 = p0 * gdScalar(11.0);

const gdScalar r0 = p0 / T0 / R;
const gdScalar r1 = p1 / T1 / R;

SodProblem::SodProblem(gdScalar gamma, gdScalar cfl, unsigned int numThreads, BOOL saveFull)
: ProblemN<1>(SODS_SHOCK_TUBE, gamma, cfl, numThreads, saveFull)
, m_xMax(0.0)
{
}

SodProblem::~SodProblem()
{

}

void SodProblem::initGrid(gdScalar xMax, unsigned int numCells,
							 gdScalar tolerance, unsigned int stepsThreshold)
{
	m_dimensions[0] = numCells;
	m_xMax = xMax;
	gdScalar dx = xMax / gdScalar(numCells);
	unsigned int numFaces = numCells - 1;
	GodunovSolverMT<1>::initModel(numCells, numFaces, dx, tolerance, stepsThreshold);
}

unsigned int SodProblem::getNumPlotTypes() const
{
	return 5;
}

std::string SodProblem::getPlotTypeName(unsigned int id) const
{
	static const std::string names[] = { "Density", "Pressure", "Temperature", "Velocity", "Momentum" };
	_ASSERT(id < 5);
	return names[id];
}

void SodProblem::plotSolution(CDC& dc, const CRect& dstRect, const CColorBar& clr, unsigned int timeStep, unsigned int plotType)
{
	const TimeStep<1> *ts = NULL;
	if (!getData(timeStep, &ts))
		return;
	_ASSERT(NULL != ts);
//	ProblemN<1>::Grid &gr = m_solution[timeStep].grid;
	const ProblemN<1>::Grid &gr = ts->grid;
	unsigned int gridSize = (unsigned int)gr.size();

	const gdScalar* start = NULL;
	unsigned int offset = sizeof(ProblemN<1>::Cell) / sizeof(gdScalar);
	gdScalar minVal, maxVal;

	switch(plotType)
	{
		case 0: // Density
			start = &gr[0].r;
			minVal = m_minValues.r;
			maxVal = m_maxValues.r;
		break;

		case 1: // Pressure
			start = &gr[0].p;
			minVal = m_minValues.p;
			maxVal = m_maxValues.p;
		break;

		case 2: // Temterature
			start = &gr[0].T;
			minVal = m_minValues.T;
			maxVal = m_maxValues.T;
		break;

		case 3: // Velocity
			start = &gr[0].v[0];
			minVal = m_minValues.v[0];
			maxVal = m_maxValues.v[0];
		break;

		case 4: // Momentum
			start = &gr[0].u[0];
			minVal = m_minValues.u[0];
			maxVal = m_maxValues.u[0];
		break;

		default:
			_ASSERT(0);
		break;
	}

	CRect tmp = dstRect;
	drawPlot(dc, tmp, clr, 0, m_minCellSize * gdScalar(m_grid.size()), minVal, maxVal);

	gdScalar val = *start;
	gdScalar clrFactor = (val - minVal) / (maxVal - minVal);
	int y = tmp.bottom + int(gdScalar(tmp.top - tmp.bottom) * clrFactor);
	dc.SelectObject(GetStockObject(DC_PEN));

	dc.SetDCPenColor(RGB(0, 0, 0));
	dc.MoveTo(tmp.left, tmp.top);
	dc.LineTo(tmp.right, tmp.top);
	dc.LineTo(tmp.right, tmp.bottom);
	dc.LineTo(tmp.left, tmp.bottom);
	dc.LineTo(tmp.left, tmp.top);

	dc.MoveTo(tmp.left, y);

	for (unsigned int i = 1; i < gr.size(); ++i)
	{
		gdScalar factor = gdScalar(i) / gdScalar(gr.size() - 1);
		val = start[i * offset];
		clrFactor = (val - minVal) / (maxVal - minVal);
		int x = tmp.left + int(gdScalar(tmp.right - tmp.left) * factor);
		int y = tmp.bottom + int(gdScalar(tmp.top - tmp.bottom) * clrFactor);
		dc.SetDCPenColor(clr.GetColor(clrFactor));
		dc.LineTo(x, y);
	}
	//const int w = 2;
	//for (int i = dstRect.left; i < dstRect.right; i += w)
	//{
	//	gdScalar factor = gdScalar(i) / gdScalar(dstRect.Width() - 1);
	//	int id = int(factor * gdScalar(gridSize - 1));
	//	gdScalar val = start[id * offset];
	//	gdScalar clrFactor = (val - minVal) / (maxVal - minVal);
	//	dc.FillSolidRect(i, dstRect.top, i + w, dstRect.bottom, clr.GetColor(clrFactor));
	//}
}

void SodProblem::generateGrid()
{
	unsigned int numCells = m_grid.size(), i;
	gdScalar l, r;
	for (i = 0; i < numCells - 1; ++i)
	{
		GodunovSolverMT<1>::Cell &c = m_grid[i];
		l = gdScalar(i) * m_minCellSize;
		r = gdScalar(i + 1) * m_minCellSize;
		//c.dv = gdScalar(4.0) / gdScalar(3.0) * gdScalar(M_PI) * r * r * r - 
		//	   gdScalar(4.0) / gdScalar(3.0) * gdScalar(M_PI) * l * l * l;
		//GodunovSolverMT<1>::Face &f = m_faces[i];
		//f.da = gdScalar(4.0) * gdScalar(M_PI) * r * r;
//		c.dv = m_minCellSize * m_minCellSize * m_minCellSize;
		c.dv = m_minCellSize;
		GodunovSolverMT<1>::Face &f = m_faces[i];
//		f.da = m_minCellSize * m_minCellSize;
		f.da = 1.0;
		f.n[0] = gdScalar(1.0);
		f.fc[0] = i;
		f.fc[1] = i + 1;
	}

	// Last cell. Doesn't have face with the next one
	GodunovSolverMT<1>::Cell &c = m_grid[i];
	c.dv = m_minCellSize;
	//c.dv = gdScalar(4.0) / gdScalar(3.0) * gdScalar(M_PI) * r * r * r - 
	//	   gdScalar(4.0) / gdScalar(3.0) * gdScalar(M_PI) * l * l * l;
}

void SodProblem::generateInitialDistributions()
{
	unsigned int numCells = m_grid.size(), i;
//	unsigned int grad = numCells;

	for (i = 0; i < numCells / 2; ++i)
	{
		GodunovSolverMT<1>::Cell &c = m_grid[i];
		c.p = 1.0f;
		c.r = 1.0f;
		c.v[0] = gdScalar(0.0);
		c.u[0] = gdScalar(0.0);
	}

	for (i = numCells / 2; i < numCells; ++i)
	{
		GodunovSolverMT<1>::Cell &c = m_grid[i];
		c.p = 0.1f;
		c.r = 0.125f;
		c.v[0] = gdScalar(0.0);
		c.u[0] = gdScalar(0.0);
	}
/*
	unsigned int grad = numCells / 4;
	for (i = numCells / 2 - grad; i < numCells / 2 + grad; ++i)
	{
		GodunovSolverMT<1>::Cell &c = m_grid[i];
//		gdScalar factor = gdScalar(i) / gdScalar(grad - 1);
		gdScalar factor = (gdScalar(i) - gdScalar(numCells) / 2.0) / gdScalar(grad);
		factor = abs(factor);
		if (factor > gdScalar(1.0))
		{
			factor = gdScalar(1.0);
		}
		if (factor < gdScalar(0.0))
		{
			factor = gdScalar(0.0);
		}
		factor = gdScalar(1.0) - factor;
		//c.p = p0;
		//c.r = r0;
		c.p = p0 + (p1 - p0) * factor;
		c.r = r0 + (r1 - r0) * factor;
		c.v[0] = gdScalar(0.0);
		c.u[0] = gdScalar(0.0);
	}
*/
	//unsigned int grad = numCells / 10;
	//for (i = 0; i < grad; ++i)
	//{
	//	GodunovSolverMT<1>::Cell &c = m_grid[i];
	//	c.p = p1;
	//	c.r = r1;
	//	c.v[0] = gdScalar(0.0);
	//	c.u[0] = gdScalar(0.0);
	//}
}

void SodProblem::applyBoundaryConditions()
{
	unsigned int numCells = m_grid.size();
	m_grid[0].r = m_grid[1].r;
	m_grid[0].p = m_grid[1].p;
	m_grid[0].u = -m_grid[1].u;

	m_grid[numCells - 1].r = m_grid[numCells - 2].r;
	m_grid[numCells - 1].p = m_grid[numCells - 2].p;
	m_grid[numCells - 1].u = -m_grid[numCells - 2].u;
	//m_grid[0].r = m_grid[1].r;
	//m_grid[0].p = m_grid[1].p;
	//m_grid[0].u = m_grid[1].u;

	//m_grid[numCells - 1].r = m_grid[numCells - 2].r;
	//m_grid[numCells - 1].p = m_grid[numCells - 2].p;
	//m_grid[numCells - 1].u = -m_grid[numCells - 2].u;
}

void SodProblem::solveRadiation(gdScalar timeDelta)
{
	// Here you should solve radiation problem; All
	// the variables in the m_grid, such as Density r,
	// Energy e, Temperature T are up-to-date and valid
	// I GUARANTEE ;)
}